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ABSTRACT 


In the area of computational acoustics, procedures which accurately pre- 
dict the far-field sound radiation are much sought after. A systematic de- 
velopment of such procedures are found in a sequence of papers by Atassi. 
et. al. [3] [4] [7]. The method presented here is an alternate approach to 
predicting far field sound based on simple layer potential theoretic methods. 
The main advantages of this method are: it requires only a simple free space 
Green’s function, it can accomidate arbitrary shapes of Kirchoff surfaces, 
and is readily extendable to three-dimensional problems. Moreover, the pro- 
cedure presented here, though tested for unsteady lifting airfoil problems, 
can easily be adapted to other areas of interest, such as jet noise radiation 
problems. Results are presented for lifting airfoil problems and comparisons 
are made with the results reported in [3] [4] [7]. Direct comparisons are also 
made for the flat plate case. 
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1 INTRODUCTION 


The prediction of far-field sound radiation is a key area of interest in the com- 
putational aeroacoustics community. This subject is a common intersection 
for most unsteady external aerodynamic problems such as the gust response 
of airfoils, flutter problems , and jet noise. Since these problems are posed in 
open domains, in general, it is difficult to extend the computational domains 
to the far-field due to the prohibitive cost of the computations involved. 
While the state of the art of computers provide high speed and large storage 
capacity, the numerical algorithms themselves are known to have problems. 
They are subject to their dissipative and dispersive properties and may con- 
taminate acoustic calculations due to the disparity in magnitude between the 
acoustic pressure and the dominant flow quantities. Current efforts involving 
the direct numerical simulation of far-field sound are being made by Tam et. 
al. [11], Roe et. al. [8], and Mitchell et. al. [6]. Among all the numerical 
development axe a sequence of works developed by Atassi and his associates 
[3] [4] [7]. A summary of these efforts are found in [1]. These methods 
employ numerical simulations in the near field, developed by Scott [9] [10], 
and construction of a Kirchoff surface for the prediction of far-field sound 
using Green’s function techniques. Thus the driving philosophy is accurate 
numerical simulation in the near-field and a “semi-analytical” approach to 
predict the far-field sound. This paper is motivated by these works. The key 
difference here is in the formulation of the far-field sound radiation calcula- 
tions. In contrast to the work of Atassi et. al. [3] [4] [7], a potential theoretic 
approach is implemented. This method requires only the free space Green’s 
function and involves an unknown simple layer density. This simple layer 
density is solved for by a technique proposed by Hariharan and MacCamy 
[5] for electromagnetic scattering problems. The application of this method 
for aerodynamic problems is shown here in detail. The method has the key 
advantage that the Kirchoff surface does not have to be a circle or other 
simple shape. As long as it is a smooth surface, the method will work. This 
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is attractive to those studying jet noise problems, where the computational 
domains are typically elongated. Also, the method is extendable to fully 
three-dimensional problems with little difficulty. 

The calculations presented here are mainly driven by the comparison of 
the results to those of Atassi [3] and Patrick [7] in an elliptical domain. We 
also compare results in a rectangular domain for a specific flat plate case pre- 
sented at the ICASE/LARC Workshop on Computational Aeroacoustics [2]. 
The near-field calculations are performed by the GUST3D code developed 
by Scott [9] [10]. 

2 FORMULATION OF THE PROBLEM 

The problem discussed here deals with the calculation of acoustic pressure 
in the exterior flow field Q, of some boundary P, where P is a Kirchoff surface 
which surrounds an arbitrary airfoil. The airfoil is submersed in a subsonic 
compressible flow field with a three-dimensional upstream vortical distur- 
bance. The disturbance has magnitude and frequency vectors given by a and 
k respectively and the details of this disturbance are outlined in reference 
[9]. We assume that the acoustic pressure on T is known. The governing 
continuity and momentum equations for compressible subsonic flow are as 
follows: 

pt -1- div(pXJ) = 0 (1) 

U t + (U. V )U+-VP = 0 (2) 

P 

and the state equations relating p and p are 

dp 

p = Ap" 1 , — = Cq (7 = 1.4 for standard air) (3) 

dp 

where p is density, p is pressure, U is the fluid velocity, t is time, A is a 
constant, and cq is the speed of sound. We assume that T is far enough from 
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the airfoil so that the mean flow quantities do not differ appreciably from 
the free stream in Cl and linearize (1) and (2) above using 

U = U ro + u 


P = poo + p' 

p = Poo + ?'• (4) 

Uoo, Poo, and poo are the freestream velocity, density, and pressure respec- 
tively, and similarily, u, p' and p' are the perturbation velocity, density and 
pressure. Before proceeding, we note that the following relation may be 
derived using (3) and the linearization in (4): 

P' = clp'. (5) 


Substituting (4) into (1) and (2), one obtains the linearized continuity and 
momentum equations: 

j^p' + poo div u = 0 (6) 

+ = 0 (7) 

where — J-j + ^oo^-- Now, we substitute (5) into (6) and (7) to obtain 


1 D 0 , , 

V J-.. P + poodlVVL = 0 

ci, Dt 


Do , 

poo Dt U + VP = °- 

Manipulating (8) and (9), namely 7^(8) —div (9), yields 


J_J1 

cL Dt 2 


p' - V V = 0 


or 


£<S + - A * 


( 8 ) 

( 9 ) 


( 10 ) 
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Now, we non-dimensionalize (10) as follows: 


xi,® 2j *3 by c/2 

t by c/2£7oo 
h,k 2 ,k 3 by 2/c 

p' by pMa\ 


where c is the airfoil chord length. Then the dimensionless form of (10) is 

M l(ft + 9“)V = vV- (11) 

We then make the following transformation which removes the depen- 
dence on time t, and the spanwise coordinate x 2 : 

p' = p(x u x 2 )e- iklt+ik3X3 e- iKlXl . 


The nondimensional wave number ki is called the reduced frequency and is 
given by ki = u;c/2C/ 00 (the dimensional k x is given by ki = u;/[/ 00 , where u; 
is the dimensional angular frequency). The quantity Ki is defined by 

Ki = kl ™°° , where /?«, = y/l - M*,. 

Poo 


(11) then becomes 


#33 + 

+ pi 


~ h \)p — o. 


( 12 ) 


Finally, we transform the coordinate system to the well-known Prandtl- 
Glauert plane using the linear transformation 


a?i = ®i 

2 2 = A »®2 

*3 = fl<x x 3- 
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Since the quantity (3^ is less than unity, this amounts to a compression of 
the ® 2 and *3 coordinates. After applying this transformation, we have 

VP+K 2 p= 0, (13) 

where 


Now the problem clearly becomes a classical exterior problem governed by the 
Helmholtz equation (13) subject to numerically calculated values of p on the 
boundary f, where T is the transformed boundary in the Prandtl-Glauert 
plane. Even though it appears to be a simple problem, the assembly of 
solutions plays a major role in the accurate prediction of far field calculations. 
A procedure that has a general structure is proposed here and is described 
next. 


2 _ 

ft ft 


3 SOLUTION PROCEDURE 

The solution procedure is based on simple layer potential theory [5]. The 
advantage of the procedure is that the solution relies on a free space Green’s 
function rather than a specific Green’s function which is suited for the domain 
of the problem. In particular, arbitrarily shaped Kirchoff surfaces may be 
used. The free space Green’s function for (13) is 

G(i,y) = -jff 0 1 (ir|x-y|), (14) 

where x is a point in the transformed domain, Q, and y is a point on T. 
From potential theory [5], the solution of (13) can be written as 

?(x) = cr(y)G(x,y)ds£ x € fi, y € f (15) 

where p satisfies 

• p ~ as |x| — ► 00 with R = |x — yj 
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( 16 ) 


• p = /(x) on f . 

Imposing the boundary condition, (15) becomes 

J„ o-(y)G(5c, y)ds s = /(5c), 5c, y € f . 

We can now describe a procedure for solving (16) for the a values. Upon cal- 
culating cr and substituting in (15), p at any point 5c £ fi can be determined. 

To allow an arbitrary shape of the Kirchoff surface, it will be assumed 
that f is some polar representable geometry. Therefore, it can be expressed 
as the function r = 12(0). Let 0 and <f> be polar representations of the points 
x and y respectively. Then 5c and y may be expressed in polar coordinates 
as 


5c(0) = (12(0) cos 0, 12(0) sin 6) 
y(<j)) = (R((f>) cos (f>,R{(f>) sin <f>), 

and we may write 

/(5c) = f{R{0) cos 0, R{0) sin 0) = f(0), 

and 

o-(y) = <t{R{4>) cos <f), R(<p) sin <f> ) = a{4>). 

The distance between the two points 5c and y is then given by 

|5c — y| = |(12(0) cos 0,R{0) sin 0) — (R{4>) cos 4>,R(<f>) sin (f)\ 

= y/mT) + 12 2 (0) - 212(0)12(<£) cos(0 - 0) 

= d(0,<f>), (17) 

and the Green’s function becomes 

<?(*.?) = -jitf'Ws-yl) 

= H^\Ki{6A)) 

= (is) 
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Now, since y= (®,y), wtere x and y are given by 

x = R(<f>) cos (f> 
y = R{<f>) sin <f > , 

dx and dy are given by 

dx = [— R(<f>) sin <f> + R'(<f>) cos <f>] d(f> 
dy = cos 4> + R'(<f>) sin (f>] d(j) 

and ds is as follows: 

ds = yfdi 2 + dy 2 

= \J[R '( <£) cos ~4> — f?( 0 ) sin 0] 2 + [i 2 '(^) sin + f£( 0 ) cos ^] 2 

= + JPM 

Substituting the above quantities into (15) yields 

/(«) = £ <rMG(6, 4>)yjR*(4>) + JP(fl #, e € [ 0 , 2 tt ). 

By defining a(0) = a(4>)\/R' 2 + f2 2 , we obtain our final equation 

/(«) = r ?(«<?(«, m, « € [ 0 , 2 ir). (19) 

«/o 

We now consider the case where 0 = <f> in (18) and G becomes singular. 
We must first note that the zero order Hankel function of the first kind may 
be expressed by 

~\ H o\ z ) = ^log(z) + 7- J + £(*), 

where 7 = ( 7 is Euler’s constant and i?(z) — > 0 as z — > 0). Using 

this, G may be written as follows: 

G(M) = ~H^\Ki{6,4,)) 

= ^iog(iri(e,^))+7- j + ji(jr<((M)). (20) 
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Our next step is to observe that when 6 —> <f>, R{6) can be expanded in a 
neighborhood of 0 by a Taylor series in the following way: 

m = rm + rw - *) + ^k> - ■ ■ • • 

Next, we square (17) and substitute the above for R(8) to obtain 
<f-(6,4>) = iR 2 ^) + 4fl(-J)i!'(?)sin ! (^-~)(« - tj>) 

+[*"(« + sin J (^)](« - 

+ lRWR"‘W + • • • . 

Now we divide both sides by the first term on the right hand side and take 
the limit as 6 — » <f> to obtain 

- .P(M) = , , 

fi — ^ AR 2 (6) sin 2 (^|^) ' R 2 {4 > ) 

Thus, when 6 — * <f>, 

d(6,4>) = 2sin(^=-^) v /S^j + #(*). (21) 

This confines the singularity which arises in (20) to the factor sin(^^-). 

Our next step is to rewrite (19) as follows: 

- *{0))G{e, <f>)d<f> + <r(d) J* G{8, <f>)d<f> = g(6), 6 € [0, 2tt). (22) 

The motivation for doing this is that when 6 — > ^ , it can be shown [5] that 
the first term in (22) will approach zero, and we may use (21) to work with 
the singularity in the last term of (22). Substituting (20) in the last term of 
(22) gives 

ms?G{e,m 

= m + R 2 W *>[ I T a l) + i - j + *(M)» 

= 9($) J*’ Iog(sin(M)]# + 9(6) R(6, i)ii, 6 S [0, 2i) (23) 
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where R is defined by 


5 = j 7 — i + h iog(2jf + &m, 6 = $ 

\ <*(M)) ~ ^ iogfsinl!^^]), 8^4 1 

Now, we must turn our attention to the first term in (23). Using complex 
analysis, the exact value of this integral may be obtained, namely 


/ log[sin(- 
J o 


-)]d(f> = — 2-7rlog(2). 


Using this result in (22) gives us the following: 

[ 2 * (*(</>) - a(6))G(6, 4>)d(f> 

J 0 


/*2tt _ 

— a(6) log(2) + a(6) j ^ R(6,<f>)d<f> = g(6), 6 6 [0, 27t). (24) 

With this equation, we have removed the problem of the singularity in (19) 
and are now ready to attempt to numerically find the a values. 


4 NUMERICAL PROCEDURE 

To numerically determine the solution, we must first choose a suitable shape 
for the Kirchoff surface. The choice of an ellipse, specifically a circle, will 
greatly simplify the above formulation since analytical representations for 
R(8 ) and its derivative are well known. If we choose some other shape, we 
may have to numerically approximate R'(8), thereby reducing accuracy. The 
Kirchoff surface must also be far enough away from the airfoil so that the 
mean flow quantities on f deviate only slightly from those of the free stream. 
We may then use the GUST3D code to determine the values of p at a finite 
number of points on the boundary. 

A numerical representation of (24) can be obtained using rectangular 
quadrature. Letting a — 0, b = 27T, hj = <pj — and h{ = 9i — 
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equation (24) becomes 


£[*(<k) - <r{8))G{6,(f>j)h j+ 1 


i=i 


-<K*)log(2) + mil WMhj+i = />). 
i=i 

Noting that both 0 and (f> correspond to points on T, we may replace 8 by 0; 
and obtain 

i= i 

- »(«<)1o 6(2) + 5(«i) £ m.Wi+l = W) (25) 

i=l 

For each i from 1 to n — 1, an equation of the above type will be yielded. We 
will have n — 1 equations and n — 1 unknowns(the a values), and therefore 
it would be to our benefit to represent (25) as 

S Ai & = & 

3 = 1 

where ,/4y = .4(0,-, <pj), aj = and /,• = With (25) in this form, 

we could use a matrix solver to solve the augmented matrix yielded by the 
system of equations. After some scrutinization of (25), it can be seen that 
the values used for Aij should be as follows: 


-log(2) -f [Efc=i Rikhk +1 - Gikhk+i], i = j 

Gijhj+i, i ^ j 


where Rij = R(8i,<f>j) and Gij = G(8i,(f>j). 

Once our goal of finding the a values is reached, it is a rather simple 
matter to obtain p at any point x € ft. Since our solution procedure was 
derived in terms of T and ft, in order to find the solution at some point x € ft 
we must first transform it to a point x G ft. Once we have a point x G ft, we 
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can use the solution procedure outlined above. For simplicity of notation, let 
(r, 8) and (r, <f>) be polar representations of the points x and y respectively, 
and proceed as follows. Let x € be represented as 

x = (r cos 6,r sin0) 

where r and 8 are fixed, and let y € f be defined as follows: 

y = (R(<j>) cos <f>, R{4 > ) sin <j>), <t> € [0, 2ir). 

Therefore, the distance between x and y is given by 

|x — y| = ^/[r cos <f) — J2(<^) cos (j)} 2 + [r sin 8 — R(4>) cos <f)] 2 
= yfr* + i2 2 (<£) - 2ri?(<£) cos {8 - J) 

and the Green’s function becomes 

As follows from (20), the solution will be 

i>M) = -7 Hr , «,«)#. 

4 Jo 

This can be expressed numerically as 

f(r , «) = - j E <Kr, 6, 

* i=i 

5 RESULTS AND DISCUSSION 

All calculations use the above described method, where the GUST3D code 
is used to predict the pressure on the Kirchoff surface. A circular boundary 
centered at the center of the airfoil is used and the pressure values inter- 
polated from those provided by GUST3D using cubic spline interpolation. 
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The number of data points on the Kirchoff surface is a function of reduced 
frequency; from 40 points for lower frequencies to 80 for higher frequencies 
(see references [9] and [10] for details). The system of equations generated is 
solved using Gaussian elimination with pivoting. The radius of the Kirchoff 
Surface is 1.5 times the semichord for flat plates. For airfoils with thickness, 
the radius is the minimum radius for which the mean flow and ffeestream 
velocities differ by less than 5% and typically varies from 1.3-2. 5 times the 
semichord. Increasing or decreasing this radius by a small amount does not 
significantly alter the solution as is indicated in figure 15. 

5.1 ICASE Benchmark Problem 6 

In this case we have a flat plate in a transverse gust (see figure 2.b) defined 
as 

, = 0.1csin [!(-£--*)] 

where the normalization for the velocity is with respect to 0 ^; length with 
respect to Aa; = 1; and time, with respect to Ax/c^. The Mach number is 
given as 0.5, and the chord length is 30 units. Using our non- dimensional- 
ization, we get ki = 157r/4 = 11.781 and o 2 = 0.2. The unsteady pressure 
is to be calculated on a box surrounding the flat plate as shown in figure 
1. The sides of the box are located at dimensional positions of a?i = ±95 
an d x 2 = ±95. When non-dimensionalized by the semichord, the values are 
x\ — ±6.333 and x 2 = ±6.333. 

The rms pressure, (||p , | 2 ), is determined on the sides of the box and 
compared to those obtained by Patrick and Atassi [2] using a semi-analytical 
approach in figures 3-6. 100 data points per side were used. The results 
compare quite well, considering that the acoustic energy is such a small 
fraction of the total flow energy, and highly accurate results are difficult to 
obtain. 
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5.2 Polar Plot Comparisons 

Polar plots of the solution on a circle of radius r were made, where a value 
of 100 times the semichord was used for r. This value is quite arbitrary and, 
in fact, the solution was multiplied by to remove the dependence on r. 
The solution was determined at 200 evenly spaced data points on the circle. 

Direct comparisons to those of Atassi and Patrick [7] were made for the 
flat plate case for reduced frequency values of 1.0 and 3.0 and Mach numbers 
0.1, 0.5 and 0.8 in a transverse gust (see figure 2.b). These are compared in 
figures 7-12 and only very slight differences are observed. 

In figures 13-14, we examine the effects of airfoil thickness for symmetric 
airfoils at no angle of attack in a transverse gust for the Mach numbers and 
reduced frequencies given above. Several examples of polar plots of the far 
field solution for airfoils with camber and angle of attack in both transverse 
and oblique (see figure 2.b) gusts are shown in figures 16-20. These serve as a 
demonstration of what kind of information can be gleaned from this method. 
As the purpose of this paper is to simply present the potential theoretic 
approach to determine fax field acoustics as an alternative to other Green’s 
functions methods, we will not go into any indepth analysis of these plots. 
However, the expected results are obtained for these cases (see references [3] 
and [7]) with only minor discrepancies, although no direct comparisons to 
the results of Atassi and Patrick are available at this time. 

6 CONCLUSION 

It was shown that potential theoretic methods are a viable alternative to 
the modified Green’s function approach of Atassi for obtaining the far field 
acoustic radiation from lifting airfoils. The only numerical techniques needed 
are rectangular quadrature, Gauss-Jordan elimination with pivoting, and 
cubic spline interpolation- all of which are simple, well-established techniques. 
Morever, the method does not require the development of a new Green’s 
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function but uses the well-known free space Green’s function. This allows 
for easy extendability to completely three-dimensional problems, i.e., the 
solution is represented as a surface integral and the three-dimensional free 
space Green’s function for the Helmholtz operator is used. Future work is 
planned along these lines. 
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